Evolutionary dynamics, intrinsic noise and cycles of co-operation 
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We use analytical techniques based on an expansion in the inverse system size to study the 
stochastic evolutionary dynamics of finite populations of players interacting in a repeated prisoner's 
dilemma game. We show that a mechanism of amplification of demographic noise can give rise to 
coherent oscillations in parameter regimes where deterministic descriptions converge to fixed points 
with complex eigenvalues. These quasi-cycles between co-operation and defection have previously 
been observed in computer simulations; here we provide a systematic and comprehensive analytical 
characterization of their properties. We are able to predict their power spectra as a function of the 
mutation rate and other model parameters, and to compare the relative magnitude of the cycles 
induced by different types of underlying microscopic dynamics. We also extend our analysis to the 
iterated prisoner's dilemma game with a win-stay lose-shift strategy, appropriate in situations where 
players are subject to errors of the trembling-hand type. 

PACS numbers: 02.50.Le, 05.10.Gg, 02.50.Ey, 87.23.Kg 



I. INTRODUCTION 



Traditionally, modellers in biology and related disci- 
plines use deterministic ordinary or partial differential 
equations to capture the quantitative behavior of dynam- 
ical systems in those fields. Such an approach is valid 
and accurate only if stochastic effects induced by exter- 
nal or intrinsic fluctuations can be neglected. External 
noise might result from environmental factors or as an 
attempt to include the effects of numerous, but weak, 
external effects. Intrinsic fluctuations arise from the dy- 
namics of the system itself. One of the most common 
sources of such stochasticity in biology is discretization 
noise in systems composed of a finite number of interact- 
ing individuals. While deterministic descriptions can be 
derived, and shown to be exact in the limit of infinite sys- 
tem size, finite systems retain an intrinsic randomness, 
sometimes referred to as demographic noise Such 
fluctuations can invalidate conclusions based on the anal- 
ysis of the deterministic dynamics, turning deterministic 
fixed points into stochastic quasi-cycles, inducing helical 
motion about limit cycles |2|], or giving rise to Turing 
patterns induced by intrinsic noise 0, . The existence 
of stochastic quasi-cycles has been known for a number 
of decades in the context of predator-prey-like systems, 
and methods have been devised to distinguish them from 
noisy limit cycles [f| ■ 

Only very recently have systematic methods, based on 
a system-size expansion of the master equation describ- 
ing the microscopic stochastic processes, been devised to 
study them analytically [6|. These methods use an ex- 
pansion in the inverse system size Q, and are now be- 
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ing applied to a number of fields in which quasi-cycles 
have been reported, including epidemiology |8l-ll0l| , bio- 
chemical reactions (Tlj . gene regulation [l2j . and more 
recently learning algorithms of interacting agents [l3| . 
The purpose of the present work is to apply these ideas 
to problems in evolutionary game theory, and to provide 
an analytical characterization of stochastic quasi-cycles 
found in computer simulations of populations of inter- 
acting players [T3 |. 

Evolutionary dynamics in this context is a mathemat- 
ical framework describing co-evolving populations. It is 
the main tool-kit used in attempts to reconcile the evolu- 
tion of co-operation with Darwinian natural selection — 
a problem which was listed as one of the 25 most press- 
ing scientific challenges in Science magazine in 2005 [l5[ . 
The problem of how mutual co-operation is sustained in 
a population subject to selection pressure favoring self- 
ish behavior is most commonly modeled using the pris- 
oner's dilemma (PD) game |iall2|. The PD is a classic 
game-theory problem in which two players have to si- 
multaneously choose whether to co-operate or to defect. 
Although the payoff for mutual co-operation is higher 
than that for mutual defection, the payoff for defecting 
when the other player co-operates is higher still. De- 
fection then forms the Nash equilibrium of the game, i.e. 
the outcome one may expect if the interacting players are 
fully rational. A number of experiments have been per- 
formed in behavioral game theory (examples are [l8lll9j) 
and biological realizations of the PD include the s tudy 
of competitive interaction among viruses, see e.g. |20l |. 
An extension of the basic PD game considers repeated 
interaction of a given pair of players. The space of avail- 
able strategies then becomes too large to allow for an 
exhaustive analysis. Most studies therefore focus on a 
selected set of strategies, such as always defect (A11D), 
always co-operate (A11C), tit-for-tat (TFT) or win-stay 
lose-shift (WSLS). A11C players always co-operate in any 
iteration, and similarly A11D players always defect. TFT 
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co-operates in the first round and then copies its op- 
ponent's previous move. This strategy emerged as the 
winner in a computer tournament run by Axelrod in 
1981 [16|. Since then TFT has been the subject of a large 
body of work [HIIlHil . Even though TFT won a subse- 
quent second competition as well, TFT is not perfect. In 
more realistic situations where players can make mistakes 
TFT can become locked into patterns of alternative co- 
operation and defection with another TFT player [24| . It 
is also vulnerable to invasion from co-operators via neu- 
tral drift. Nowak and Sigmund [25| then proposed WSLS; 
this strategy has none of the above disadvantages. WSLS 
co-operates in the first round and then keeps playing the 
same action (co-operate or defect) if it receives a favor- 
able payoff, and switches from one action to the other 
if it does not. It can resist neutral drift by co-operators 
and can correct mistakes, avoiding disadvantageous cy- 
cles. There is evidence to suggest that some animals 
employ these strategies, for example, in their behavior in 
the presence of predators [H, [2?} ■ 

Historically, the analysis of evolutionary dynamics has 
mostly been based on deterministic replicator dynam- 
ics |28|, explicitly excluding stochastic effects. More re- 
cently, methods from statistical physics and the theory 
of stochastic processes have been used to study games 
in finite populations. In the absence of mutation, a fi- 
nite population will always fix on a given strategy due to 
stochastic fluctuations. The resulting fixation probabili- 
ties and average fixation times can be calculated [29l - l3ll | . 
Further quantities of interest are stationary distributions 
of the underlying stochastic processes [H, l32l - [34j ]. and 
the phenomenon of dynamic drift (35j . 

In the context of these studies of stochastic processes in 
game theory, cyclic behavior has been reported [36l - l38l | in 
the rock-papers-scissors game, and in [14| , where stochas- 
tic quasi-cycles between co-operation and defection have 
been observed in finite populations of agents playing the 
iterated PD. In the present work we will focus on the 
latter game, and provide an analytical theory which al- 
lows one to compute properties such as power spectra, 
or equivalently the correlation functions of these quasi- 
cycles, to a good approximation in the limit of large, 
but finite populations. Based on this analytical approach 
we are able to identify regions in parameter space where 
stochastic quasi-cycles would be expected to occur, and 
we compare the amplitude of cycles arising from different 
types of microscopic update dynamics. 

The paper is organized as follows. In Sec. [TT| we out- 
line the iterated PD and define the different microscopic 
processes. We first focus on the case of only three pure 
strategies, A11C, A11D, and TFT. The deterministic anal- 
ysis for this model is presented in Sec. Mil with a classi- 
fication of the fixed points and an exploration of the pa- 
rameter space. We move from a deterministic description 
to a stochastic formulation in Sec. lIVI and consider effects 
arising in finite populations. In particular, we carry out 
a system-size expansion of the master equation allowing 
us to classify the periodic stochastic deviations from the 



deterministic limit. In Sec. |V]we extend the analysis to 
include WSLS as a fourth strategy. Finally, in Sec. I VII 
we summarize our findings and outline avenues of future 
research. 



II. MODEL AND DEFINITIONS 
A. Iterated PD 

We will mostly follow the setup of Imhof et al 14]. An 
exception will be when we discuss the extension of the 
model in Sec. [V] As such, we will consider a popula- 
tion of N players with each player carrying out one of 
three pure strategies: A11C, A11D, or TFT. The respec- 
tive payoffs resulting from an encounter of two players is 
characterized by the following payoff matrix: 

AUC AUD TFT 

AUC ( Rm Sm Rm \ 

AUD j Tm Pm T + P(m-1)\, 

TFT \Rm-c S + P(m - 1) - c Rm - c J 

(1) 

where m is the number of rounds played when two play- 
ers meet. The parameters T,R,P, and S are the payoffs 
of the basic PD game (in which players meet only once) : 
T is the temptation to defect, i.e. the payoff a defec- 
tor receives when playing a co-operator, R is the reward 
for mutual cooperation, P is the punishment for mutual 
defection, and S is the sucker's payoff for co-operating 
with a defector. The so-called complexity cost, c, is im- 
posed on the TFT strategy and represents the alloca- 
tion of resources used to remember an opponent's last 
move [Til [39l ] . For the dilemma to be present we require 
that the parameters satisfy T > R > P > S and also that 
R > (T + S)/2, to prevent mutual alternate co-operation 
and defection being more profitable that of mutual coop- 
eration [l6|. Throughout this paper we use the specific 
parameter values T = 5, R = 3, P = 1, S = 0.1, and 
m = 10 [lij]. In the terminology of game theory, the 
iterated PD as defined by the above payoff matrix is a 
non-cooperative symmetric game. 

In the following we will label the strategies A11C, AUD, 
and TFT by i = 1, 2, 3 respectively. The number of play- 
ers in the population using strategy i will be denoted by 
rii, and we require that ri\ + + = N. The expected 
payoff, or fitness, of a player of type i is then given by 

*< = N-l ' ( 2 ) 

where the elements of the payoff matrix, Eq. (fTJ , 

e.g. an = Rm, an = Sm, etc. In using the defini- 
tion ^ we follow the choices of [TtJ and exclude inter- 
actions of one individual with itself. Definitions with 
self-interaction are possible, the differences do not affect 
the results to the order in inverse system size we will be 
working at. 
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The so-called reproductive fitness of an agent carrying 
pure strategy i, /j, is denned as [17[ 

fi = 1 -W + WTTi, (3) 

where w is a selection strength that determines the im- 
pact that the game has on the agent's overall fitness. If 
w = 0, then fi — 1 for all i, and one recovers the limit of 
neutral selection. For w > selection becomes increas- 
ingly frequency dependent. The average reproductive fit- 
ness in the population is then given by 

i 

and the average payoff is tt = J2i( n i/^) n i- in order to 
complete the model we need to specify the microscopic 
dynamics of the system, i.e. we need to define the rules 
by which the composition of the population changes over 
time. There are several such microscopic processes which 
have been studied in the literature, and we will define 
some of these in Sec. lII Cl Before we do so, it will however 
be helpful to discuss the standard replicator-mutator dy- 
namics commonly considered in the literature. Provided 
the microscopic dynamics are chosen appropriately, these 
equations are a suitable description in the deterministic 
limit, valid for infinite populations. It is however im- 
portant to stress that the replicator-mutator dynamics 
are not the limiting deterministic dynamics for all micro- 
scopic processes, as pointed out in [32],|40[, and as we will 
discuss in more detail below. 



B. Canonical replicator-mutator equation 

Within the standard replicator dynamics of evolution- 
ary game theory the time evolution of the concentration 
of a strategy i is given by [28| 

x i = Xi (fr-<t>°°), (5) 

where Xi denotes the concentration of strategy i in the 
population in the deterministic limit: Xi = lim rii/N. 

N—too 

Similarly, we will write 

f°° = lim /, = l- aj + i«y OijXj , (6) 

j 

and 

0°° = lim 4> = Y J Xifr> (7) 

N— ¥oo * — ' J 

j 

where the superscripts indicate that Eq. ([5]) are, for suit- 
ably chosen microscopic dynamics, valid only in the de- 
terministic limit of infinite populations. The basic as- 
sumption underlying these dynamics is that individuals 
reproduce asexually in proportion to their reproductive 
fitness, and that offspring inherit the strategy of their 
parent. 



If one introduces mutation, so that there is a finite 
chance that a player will produce an offspring which does 
not use the same strategy as their parent, the above dy- 
namics needs to be modified, and the description is then 
in terms of so-called replicator-mutator equations [4l| . 
Focusing on the case of M pure strategies we will as- 
sume that in a reproduction event a player produces an 
exact copy of itself with probability I — (M — l)u and 
a mutant which plays one of the other M — 1 strategies, 
each with probability u. The parameter u is confined to 
the physically meaningful range < u < 1/M for the 
case of M pure strategies. For u — 1/M an offspring will 
be of any of the M types with equal probability 1/M. It 
is convenient to introduce a mutation matrix 

The replicator-mutator equation is then given by [4l| 

Y, r '-C<h, (Q) 

In the limit of zero mutation Eq. ([9]) reduces to the stan- 
dard replicator equation, Eq. ([3]). 

C. Microscopic dynamics 

We will now define the different microscopic processes 
we will consider. We will restrict ourselves to dynamics 
conserving the total number of players in the population. 
To specify a process it is then sufficient to define the 
'conversion' rates T^j, corresponding to events in which 
a player of type i is replaced by one of type j. For the 
general case with M pure strategies i, j — I, . . . , M. We 
will limit the discussion to processes of the general form 

T i-+3 = X! Jf-j$9ki(f)qkj, (10) 
k 

where / = . • • , /m)- The form (fTU|) is found by, 
at each time step, selecting two players, one for poten- 
tial reproduction and one for potential removal, from the 
population. The player selected for potential removal is 
assumed to be of type i, and each term in the sum cor- 
responds to selecting a player of type k for potential re- 
production. A given combination (i, k) thus occurs with 
probability (riink)/N 2 . Here we use sampling with re- 
placement. Dynamics without replacement of an already 
chosen player are equally possible, leading to, for exam- 
ple, factors of N(N — 1) in the denominator instead of 
iV 2 . The differences amount to effects of order TV -1 , and 
do not affect results to the order in the inverse system 
size we will be working at. For a given pair of selected 
players reproduction and death actually only occur at a 
rate proportional to gki(f), which here we assume to be 
a function of the reproductive fitnesses (implying a possi- 
ble dependence on the average fitness 4>). The factor q^j 
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accounts for potential mutation events. The four kinds 
of microscopic dynamics we will consider in the following 
differ in the details of the function g, which we will de- 
scribe below. Choices in which g^i does not depend on / 
correspond to neutral selection. 

Before we define the details of the different microscopic 
dynamics some general statements are appropriate. For 
simplicity, we will focus on the case of a game with M = 3 
pure strategies and in particular the iterated PD game 
with strategies A11C, A11D and TFT as introduced above; 
generalization to an arbitrary number of strategies M 
is however straightforward. For any choice of gki, the 
state of the TV-player population is defined by the num- 
ber of individuals using the A11C and A11D strategies: 
n = (711,712), the number of TFT players is then given 
by 77,3 = N — ni — 7i2 • Furthermore the reproductive fit- 
nesses / and the average reproductive fitness, <fi, are fully 
determined by the state n of the system (see Eqs. (J2J)- 
(j4])). ft follows that the transition rates T^j can be 
written as functions of n, and the microscopic stochastic 
process is described by the following master equation for 
the probability, P(n,t), of the system being in state n: 



dP(n,t) 
dt 



= (Uj-i^WPKi) (11) 

+ (E 2 - l)T 2 ^ 3 (n)P(n,t) 
+ (£ 1 £ 2 - 1 -l)T 1 ^ 2 (n)P(n,t) 
+ (E 2 E^ -l)r 2 ^i(n)P(n,t) 

+ (% 1 -l)T 3 ^i(Tl)P(Tl,i) 

+ (E^-l)T 3 ^ 2 (n)P(n,t). 

Here we have introduced shift operators Ei, where i = 
1,2, acting on functions of the state of the system, 
if)(ni,n 2 ), as follows: 



Eiip(ni,n 2 ) = ip(ni + l,n 2 ), 
E^ 1 ip(ni,n 2 ) = ip(ni - l,ra 2 ). 



(12) 



Similar definitions apply for E 2 and E 2 . Multiplying 
both sides of Eq. (|11[) by n and summing over all possible 
states and using a decoupling approximation, valid for 
TV — s- 00, we can then write the rate of change of the 
concentration of strategy i as 



Xi 



[Tk^i(x)-Ti^ k (x)], 



(13) 



after a re-scaling of time by a factor N. Clearly this equa- 
tion is not restricted to the case M = 3, and holds when 
an arbitrary number of strategies are present. Transition 
rates in Eq. (Til?]) are found from Eq. (ITU|) using the sub- 
stitution 77; /N — > Xi and fa —> f°°, see Appendix A for 
further details. Substituting these limits into Eq. (fl~3)) 
one finds that the deterministic evolution of the concen- 
trations of strategies is given by 



^ ^ Xj [xkfJjkQji XigjiQjk] 



(14) 



For different update rules this equation differs only in the 
specific form of g used. 

We will now proceed to give the specific form of the 
function gki(f) for a set of different update rules which 
have previously been proposed: the Moran process, a 
linear Moran process, a local process and the Fermi pro- 
cess 132, 351. 



1. The Moran process 

In the Moran process [42I ]. once a player of type k has 
been chosen for potential reproduction and a player of 
type i for potential removal, the reproduction event oc- 
curs at a rate proportional to fk/4>i specifically we will 
choose 



h 
2<j)' 



(15) 



The arbitrary pre-factor of 1/2, equivalent to choosing 
a time scale, has been introduced to allow better com- 
parison with other update rules [35!]. By substituting 
Eq. (JTSJl into Eq. (TJU) and using Eq. Q, J2 k x k = 1. and 
J2k1j k = 1' one nn ds specifically for the Moran process 
that 



Xi 



x jfj°1ji 



Xi4>° 



2<t>° 



(16) 



It is important to stress that the average reproductive fit- 
ness <f>°° is a function of the concentration vector x, and 
so is a time-dependent quantity. While Eq. (fT6)) is 
similar to the standard replicator-mutator dynamics, the 
pre-factor (2(f>° )~ 1 corresponds to a dynamic re-scaling 
of time, and so may affect the transient dynamics. The 
location of fixed points and their local stability are how- 
ever not affected, as a straightforward calculation shows. 



2. The linear Moran process 

The linear Moran process is defined by the following 
choice [35( 



(17) 



where A > is a constant parameter, such that it is 
always the case that Tj_j.j > 0. Notice also that one 
has gk t M (f) = (1 + Att^tt/j — tt))/2. A common choice, 
which we will adopt in the following, is A = 1/Af max , 
where Af max is the maximum possible difference between 
fi and (j) [lil, i.e. Af max = max fci „ |/ fc (n) - 0(n)|. In 
the absence of mutation (u — 0) the deterministic limit 
results in the following dynamics 



2A/ m 



(18) 
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Therefore, up to a re-scaling of time by the constant fac- 
tor (2Af max )~ 1 , the linear Moran process without mu- 
tation is described by the standard replicator dynamics 
in the limit of infinite population size. However for u^O 
one does not recover the standard replicator-mutator 
equations, Eq. from the linear Moran process. 

In both Eq. (TT51) and (for u^O) from the result of sub- 
stituting Eq. (|T7)l into Eq. (ITi|) . the reproductive fitness 
only enters in differences of the form — <f> or //. — f t 
and is normalized by Af max . Since both, fitness differ- 
ences and Af max , scale linearly in w, the deterministic 
dynamics is independent of the selection strength w for 
the linear Moran process. Finally, the linear Moran pro- 
cess can be obtained from Moran process Eq. (fT5)) in the 
weak selection limit, w<l. Using Eq. ([3]) one has 



by |32j,|4S 



9ki(f) = 



1 



f 



tanh I - (/ fc - fi) 



(21) 



Unlike the other update rules, the Fermi process does not 
reproduce either the standard replicator or replicator- 
mutator equations in the deterministic limit. From 
Eq. (|2"Tj) one has 

9kiif) = ~ (l + - n)) + 0(w 2 ), (22) 

in the limit of weak selection. This is of the same func- 
tional form as Eq. (|2U|) . 



1 — W + WTTk 

2(1 — w + wtt) 

l(l + w(TT k -n))+0(w 2 ), (19) 



so that to linear order one recovers Eqs. (fl7|) with the 
choice A = 1. 



3. The local process 



III. RESULTS OF THE DETERMINISTIC 
ANALYSIS 

As an initial step towards characterizing the outcome 
of the iterated PD, we compute the fixed-point structure 
in the deterministic limit of the four different dynamics 
defined above, as function of the complexity cost c and 
the mutation rate u. We fix the selection strength to 
w = 1 throughout. 



The so-called local process was first proposed by 
Traulsen et al [iol ] and is based on a pairwise comparison 
of one agent's fitness with another in order to determine 
whether or not reproduction occurs. This process has 
the advantage that no knowledge or computation of the 
average fitness of the population is required to carry out 
a microscopic step. The local process is defined by 

where Af max is again required for normalization and 
fixed at the beginning, and then remains unchanged as 
the dynamics proceeds. As opposed to the case of the 
linear Moran process, Af max is now the maximum possi- 
ble absolute difference between any two fitnesses fi and 
f k : Af max = maxi t k,n |/*(n) - /k(n)|. As with the lin- 
ear Moran process the local process does, up to a con- 
stant factor which can be absorbed in the definition of 
time, reproduce the standard replicator equation in 
the deterministic limit if mutation is absent [HI, EOj- At 
finite mutation rates one does not however recover the 
replicator-mutator equation, Eq. ([9]) 32] . 



4- The Fermi process 

Finally, the so-called Fermi process is an alternative 
pairwise comparison process which uses the Fermi-Dirac 
distribution instead of the linear functional dependence 
on fitness differences as in the local process. It is defined 



A. General fixed point structure 

The qualitative picture one finds is similar for any of 
the four dynamics; two different threshold values of the 
mutation rate can be identified, we will refer to these 
as Uc and u c 2 K For < u < Uc\ one typically finds 
three fixed points: a locally stable attractor near A11D, a 
saddle point also near A11D, and an unstable fixed point, 
located close to the A11D/TFT edg e of the strategy sim- 
plex, see Fig. HJa). Following [1J| we will refer to this 

latter fixed point as the 'mixed fixed point'. At u = u c 
the mixed fixed point becomes stable, as shown in panel 
(b) of Fig. [TJ At u = Uc , the two fixed points near A11D 
annihilate, leaving the mixed fixed point as the only at- 

tractor for u > u c , see Fig.[IJc). At u = 1/3 the mixed 
fixed point is at or near the center of the simplex, see 
Fig. Hid). At this maximal physical meaningful value 
of u an individual of any type produces an offspring of 
any of the three different strategies with equal probabil- 
ity. While this qualitative picture is the same for all four 
dynamics considered here, the numerical values of u c ^ 

(2) 

and u c will in general be different for the different dy- 
namics, and they may also depend on the choice of the 
complexity cost, c. The overall picture is consistent with 
the results of [14[ , where the standard replicator dynam- 
ics were studied and where similar qualitative behavior 
was found. Our analysis thus demonstrates that the find- 
ings of i jyj generalize to a broader class of dynamics. The 
only difference between our findings compared to those of 
earlier analyses, lies in the saddle point described above, 
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(a) u = 10" 
A11C 




(b) u = 0.004 



AI1C 



. i i r . 
. \ i » r / . 



A11D 



TFT 



AUD 



TFT 



(c) it = O.007 



(d) u = 1/3 



FIG. 1: (Color online) Fixed-point structure and flow fields 
of the standard replicator-mutator equations for the iterated 
PD at c = 0.8 and w = 1. Black symbols are stable fixed 
points, white symbols are unstable and gray symbols denote 
saddle points. The A11D fixed point and saddle point can 
be seen in the bottom left-hand corner of the simplices. The 



(i) 



i 0.0016. 



mixed fixed point (triangle) changes stability at u. 
The AUD fixed point (circle) and the saddle point (square) 
annihilate at u; 'fa 0.005. Arrows indicate the direction of 
the deterministic flow in the strategy simplex. The color map 
shows the Euclidean speed of the trajectories, These 
images were produced using a modified version of the Dynamo 
package 



which was not reported in [Tij . presumably because it 
is not an attractor of the dynamics. Although for com- 
pleteness we have given a general account of the fixed 
point structure, the mixed fixed point will be the focus 
of our analysis in the following sections, as it is this fixed 
point which gives rise to coherent oscillations induced by 
demographic noise. 



B. Limit of small mutation rates 

Further analytical progress is possible in the limit of 
small mutation rates, u< 1. For all four cases considered 
here the deterministic dynamics, derived from Eqs. (|13[) . 
are of the form 



x = h ( o)(x) + uh (1) (x), 



(23) 



with the mutation rate entering linearly in the resulting 
differential equations. We will now make the following 
ansatz for a fixed point x*\ 



'(0) 



(24) 



< 



1< 
0.8 

0.6 

n 

I 

0.4 
0.2 




A11D 



Saddle 



Mixed 



10 



10 



10 

LI 



10" 



1 k 

(a) \ (b) 

0.98- 

0.96- 
0.94- 

o 

0.92- 

0.9- 





10 



0.005 



x 



A11C 



FIG. 2: (Color online) Panel (a) shows the xaud components 
of the three fixed points of the Moran dynamics at c = 0.8 
with changing mutation rate u. Circles are the results from 
numerically solving the fixed-point equations (|25|) . lines are 
from Eq. (|24]). Panel (b) shows the paths of the A11D fixed 
point and saddle point in phase space as u is varied. The 
fixed points meet and annihilate at u « 0.005, for u larger 
than this value neither fixed point is present. 



Here is a 



in the limit of small mutation rates u 
fixed point of Eq. (TP3"]) at u = and x%s captures the 
effect of non-zero mutation rates at next-to- leading order. 
Inserting Eq. (|24p into the fixed-point condition 



h (0) {x*) + uh (1) (x*) = 0, 
and collecting terms in linear in u, one then finds 



-J 1 h (1) {x 



(0); 



(25) 



(26) 



no- 



where J is the Jacobian of h( \ evaluated at x* 
Since a;* ^ can be found in closed form from Eqs. (TP3")) 
with u = 0, substituting, Eq. (|2"o) into Eq. (j2~4")l gives an 
analytical prediction of the location of the fixed point at 
small u. 

In Fig. [2] we compare the outcome of the above linear 
expansion with results from a direct numerical evaluation 
of the fixed points of Eq. (Tl~3"|) obtained using a Newton- 
Raphson procedure. The expansion is seen to be a good 
approximation for the location of the fixed point for val- 
ues of u up to u w 0.01. From the figure we see that the 

(2) 

A11D and saddle fixed points annihilate at u c ~ 0.005 
for this value of c. This annihilation is consistent with the 
disappearance of the A11D fixed point at large u reported 
by Imhof et al 14 [ . 



C. Mixed fixed point and phase diagram 

For suitable choices of the model parameters c and u, 
the mixed fixed point can be a stable attractor with com- 
plex eigenvalues of the corresponding Jacobian. One can 
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FIG. 3: A phase diagram showing the regions in the (c, tt) 
plane for which the mixed fixed point of Eq. (| 13f) for the 
Moran process is a stable node, stable spiral, unstable spiral 
or is orbited by a limit cycle. The solid line at u = 1/3 
indicates that at this value the mixed fixed point becomes a 
stable node situated at x — (1/3, 1/3). Note that due to the 
procedure used to check for the presence of limit cycles, the 
border between the limit cycle and unstable spiral regions is 
only approximate. 

thus expect coherent stochastic oscillations to arise in 
finite populations at those model parameters. We there- 
fore focus our attention on the mixed fixed point, and 
identify the regions in the (c, w)-plane where such com- 
plex eigenvalues are found. More generally we will deter- 
mine the nature of the mixed fixed point as a function 
of u and c. The result of numerically solving for fixed 
points of the deterministic dynamics corresponding to 
the Moran process is shown in Fig. [3l We will denote 
fixed points with purely real eigenvalues as 'nodes' and 
those with complex eigenvalues as 'spirals'. At low values 
of c we observe a re-entry phenomenon, where the mixed 
fixed point goes from a stable spiral to a stable node and 
back to a stable spiral as u is decreased. 

Numerically we also observe a region where the dy- 
namics converges onto a limit cycle. We are at this point 
unable to provide a proof for the existence of limit cycles 
or to analytically determine the position of the border 
between the limit cycle and unstable spiral regions. We 
therefore determine the presence of limit cycles by nu- 
merically integrating the deterministic dynamics using 
an Euler forward method, starting from the center of the 
simplex, allowing for a period of equilibration and then 
applying a suitable threshold criterion to detect closed 
trajectories. The unstable spiral region is identified as 
the region where we do not find limit cycles numerically. 
In situations where there is more than one attractor (e.g. 
a limit cycle and a stable fixed point near A11D) initial 
conditions will determine the stationary state of the dy- 
namics. At u = 1/3 the mixed fixed point is located in 
the center of the strategy simplex, and becomes a stable 
node. 

All other update rules studied in this paper have the 



Stable nodes 




c 



FIG. 4: Phase diagram for the Fermi process. We again see 
the same qualitative structure as for the other update rules, 
with a difference only in the positions of the boundaries. 



same qualitative features as the Moran process, and 
hence their phase diagrams are structurally similar to 
that shown in Fig. |3l except that the mixed fixed point 
does not become a stable node at x = (1/3,1/3) at 
u = 1/3 for rules that use pairwise comparison. Instead, 
the fixed point forms a stable spiral close to the center of 
the simplex. Although qualitative features of the phase 
diagrams are the same for all four update rules, the quan- 
titative positions of the borders in the (c, u) plane may 
differ for each update rule. For example, Fig.0]shows the 
stability map for the Fermi process. Here the re-entry re- 
gion persists for larger values of c and the region in which 
the mixed fixed point is unstable is also much larger. 



IV. STOCHASTIC EFFECTS AND 
SYSTEM-SIZE EXPANSION 

Until now we have focused on the dynamics of infi- 
nite populations. In this section we investigate effects 
arising in finite populations, especially stochastic oscil- 
lations arising via a coherent amplification of intrinsic 
fluctuations. Such oscillations have been found in a va- 
riety of systems as described in the introduction. These 
quasi-cycles typically arise in regions of the phase dia- 
gram where the deterministic dynamics approach a fixed 
point, and so the range of parameters in which systems of 
finite populations display oscillations is generally wider 
than the region in which the deterministic system allows 
for periodic solutions. Fig. [5] indeed confirms that this is 
also the case for the evolutionary dynamics of the iterated 
PD. In the figure we choose model parameters such that 
none of the four deterministic dynamics approach peri- 
odic attractors, but instead have stable fixed points with 
complex eigenvalues (stable spirals) . As seen in the figure 
the dynamics in finite populations still generate oscilla- 
tory behavior, induced by intrinsic fluctuations. This os- 
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FIG. 5: (Color online) The results for the concentration of the 
A11C strategy from one run of a Gillespie simulation [45|, [4(| 
for each of the four update rules at N — 10000, c = 0.8, and 
u = 0.05. The time averaged concentration of each run has 
been subtracted from the data to give the deviation from the 
deterministic fixed point. 



cillatory behavior is similar to that reported in [14j . The 
four panels demonstrate that the quality and frequency of 
these stochastic oscillations can vary over a wide range 
depending on the details of the microscopic dynamics, 
and so we will now go on to characterize their properties 
in more detail in order to obtain a more comprehensive 
picture of this phenomenon. 

The analytical approach we will use to characterize 
these fluctuations is based on an expansion of the mas- 
ter equation in the inverse system size [7]. This method 
is now standard in the analysis of interacting-agent sys- 
tems, and we will therefore not present the full details of 
the calculation, but instead restrict ourselves to giving a 
few of the intermediate steps and the final results. The 
starting point of the system-size expansion is an ansatz 
of the type 



N 



(27) 



amounting to a separation of deterministic and stochas- 
tic contributions to the number, rij, of individuals of type 
i in the population. The first term on the right, xi(t), is 
the deterministic trajectory, and captures fluctua- 
tions about this trajectory; the magnitude of these fluc- 
tuations is expected to be of order iV -1 / 2 , as reflected by 
the above ansatz. One proceeds by inserting this ansatz 
into the master equation and focuses on the prob- 
ability distribution of £, rather than that of n, so that 
one sets P(n, t) = IT(£, t). Expanding the resulting mas- 
ter equation for II(£,i) in powers of TV -1 / 2 one recov- 
ers, at leading order, the generalized replicator-mutator 
equation, Eq. f|13|) . At next-to-leading order in N^ 1 / 2 a 



Fokker-Planck equation of the form 
dU v-^ d 1 n 

i i,3 



d 2 n 



(28) 



is found 7], where Cj = J2k Jikik- Here J is the Jacobian 
of Eq. (fT3f and B is a symmetric, 2x2 matrix, whose 
precise form will depend on the exact nature of the mi- 
croscopic dynamics, but whose general form is given in 
Appendix O The Fokker-Planck equation (j2"5)l is equiv- 
alent to the linear Langevin equation 47] 

i = JZ + V, (29) 

where rj is Gaussian white noise with correlations 

(Vi{t)Vi(t'))=B ij 5(t-t'). (30) 

In contrast to the Langevin equations derived using the 
Kramers-Moyal expansion |32| . Eq. ([2U]) contains addi- 
tive, rather than multiplicative noise. In the application 
we are considering here, we are interested in fluctuations 
about the stationary state and so the matrices J and 
B are evaluated at the fixed point of the deterministic 
dynamics. 

Given the linearity of Eq. (f2U| , it is straightforward to 
compute the power spectra of the fluctuations £ about 
the deterministic fixed point. Following the steps of [llj . 
one obtains 



^ ki i 



(31) 



where $ = iuA — J and I is the 2x2 identity matrix. 

In Fig. [6] we give an example of the power spectra of 
fluctuations about the deterministic fixed point obtained 
for the Moran update rule at c = 0.8 (w = 1), u = 0.01, 
and N — 12000. Theoretical predictions from the van 
Kampen expansion and numerical simulations are in near 
perfect agreement. The spectrum shows a pronounced 
peak at a frequency of approximately uj = 0.05, indicat- 
ing the existence of amplified oscillations with that char- 
acteristic frequency. The amplitude of these oscillations 
is proportional to A^ 1 / 2 , see Eq. (|2"7|) . the proportional- 
ity constant is determined by the area under the power 
spectrum. Depending on the choice of parameter values 
one can then expect the amplitude of the quasi-cycle will 
be of order one up to system sizes of 10 or so, i.e. com- 
parable to the species concentrations at the deterministic 
fixed point. Even for very large populations the oscilla- 
tions can therefore be significant. If the trajectory of 
the system is monitored over a time scale much smaller 
than the oscillation period, then this may lead to inter- 
vals in time in which the concentration of A11C is found 
to be consistently higher than that of TFT or A11D, that 
is, to intermediate periods where co-operation dominates 
the population. Such effects may for example be relevant 
when evolutionary time scales are much longer than time 
windows over which measurements can be made. 
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FIG. 6: (Color online) The power spectra for oscillations in 
the concentrations of the three strategies with Moran updat- 
ing for c = 0.8 and u = 0.01. Symbols are the results of a 
Gillespie simulation with N = 12000 and approximately 10 4 
runs. Solid lines are theoretical predictions obtained from 
Eq. (|31l) . Simulation results show excellent agreement with 
the theory. 
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FIG. 7: (Color online) A comparison of the power spectra for 
oscillations in the A11C strategy concentration at c = 0.8 and 
u = 0.05. N ranges from 10 4 to 10 6 and the number of runs 
for each simulation is of order 10 4 . 



The magnitude of the peak in the power spectra is 
a good proxy for the amplitude of the stochastic quasi- 
cycles, and the height of the peak is in turn largely deter- 
mined by the inverse of the real part of the relevant eigen- 
value of the deterministic dynamics at the fixed point. In 
the deterministic system, perturbations about the fixed 
point decay with a time constant proportional to the in- 
verse of this real part, and it is intuitively easy to see 
that the magnitude of stochastic oscillations diverges as 
the real part of the largest eigenvalue tends to zero. More 
specifically, as shown in [48|, the magnitude of the peak 
in the spectra diverges as the system approaches a Hopf 
bifurcation, where the stable spiral becomes an unsta- 
ble one. The resulting delta-function peak in the power 
spectrum indicates that a limit cycle is born in the un- 
stable phase. This can also be seen from Eq. (piTj) and the 
definition of the matrix $. At the Hopf bifurcation the 
relevant eigenvalue of J is purely imaginary, and when u 
becomes equal to the imaginary part of this eigenvalue, 
the matrix $ becomes singular, such that the expression 
on the right-hand side of Eq. pip , involving the inverse 
of diverges. 

In order to compare the relative magnitude of stochas- 
tic oscillations in the four different dynamics at fixed val- 
ues of u and c, it is therefore useful to determine how far 
or near to the instability the pair (u, c) places the re- 
spective dynamics. In Fig. [3] we plot the instability lines 
indicating the occurrence of a Hopf bifurcation in the 
(u, c) plane for the four different types of dynamics. For 
any fixed c one finds that u[}p > u^ L > u^ LM > u^ M 
and that accordingly for any u sufficiently large to place 
all four dynamics in the stable regime, the Fermi process 
is much closer to the limit-cycle regime than the other 
types of dynamics, and would therefore be expected to 
have a larger peak in the power spectra. As discussed 
above we furthermore find that the Fermi process, with 
its alternative form of the pairwise comparison process, 
produces demographic oscillations of a higher frequency 
than the other update rules, see Fig. [7J 



V. ITERATED PRISONER'S DILEMMA WITH 
ERRORS 



Having shown that the analytical approach captures 
the properties of quasi-cycles accurately, we can now 
compare the magnitude of the stochastic oscillations for 
the different update processes at the same values of c 
and u. The power spectra of the fluctuations in the A11C 
concentration are shown in Fig. [7J for the four update 
rules at one fixed mutation rate and for a specific choice 
of the complexity cost. Results indicate that the Fermi 
process produces demographic oscillations of a higher fre- 
quency than the other update rules, in line with the time 
series shown in Fig. [5] Even though the power spectra 
for the Moran and linear Moran update rules are seem- 
ingly indistinguishable in Fig.[JJ they are not analytically 
equivalent. 



In this section we study an extended version of the 
iterated PD game, allowing for a fourth pure strategy, 
win-stay lose-shift (WSLS). It is appropriate to include 
this strategy in the discussion when so-called 'trembling- 
hand' errors are considered [23|. Trembling-hand errors 
introduce the possibility of a player making a mistake 
after they have decided what to play, that is, a player co- 
operating when they meant to defect, or defecting when 
the intention was to co-operate. We will assume that the 
two players make errors of this type independently with a 
small probability e > in any given round. TFT's disad- 
vantage is then that it can become locked into a cycle of 
alternate co-operation and defection with another TFT 
player after a mistake occurs. In such games, TFT can be 
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FIG. 8: (Color online) The borders between the stable spi- 
ral region (above the respective lines) and limit cycle region 
(below) for the different update processes. The black dot in- 
dicates the point (c, u) — (0.8, 0.05) used in Fig. 7. 



outperformed by WSLS [l7| • WSLS co-operates initially 
and then keeps using its strategy (co-operation or defec- 
tion) whenever it receives payoff T or R and switches its 
strategy (from co-operation to defection or vice versa) if 
it receives P or S. WSLS does not become locked in such 
cycles when playing against TFT or WSLS. We include 
the WSLS strategy in our game, extending the dynamics 
to three degrees of freedom. 

In the presence of trembling-hand errors the outcome 
of a PD game between two fixed players and iterated 
for a finite number of rounds will generally be stochas- 
tic and depend on the timing at which errors occur in 
the interaction sequence. In order to simplify matters 
we will therefore follow |23i] and restrict the discussion 
to cases in which an interaction between two players 
consists of an infinite number of iterations of the PD 
game. It is then appropriate to use the expected pay- 
offs per round, i.e. for two fixed players, say of types 
i,j £ {AUC, A11D, TFT, WSLS}, one formally consid- 
ers an infinite sequence of PD interactions, and uses the 
mean payoff per round to define the payoff matrix ele- 
ments ay . The payoff matrix can then be worked out for 
small error rates, and is given in [23j and reproduced in 
Appendix B for convenience. The complexity cost, c, is 
no longer a relevant parameter now that the number of 
rounds is infinite. 

Previous work on this game has shown that in the lim- 
its of zero mutation and weak selection the population 
can either fix on WSLS or A11D depending on the values 
used in the payoff matrix [23[. We continue to use the 
parameter values given in Sec. Ill Al and explore how the 
dynamics of the game depend on mutation and error rates 
and identify and classify demographic oscillations. Ana- 
lyzing the four update rules given in Sec. Ill CI we again 
find a mixed fixed point on which we focus our analysis 
— since demographic oscillations may occur about this 
fixed point when it is a stable spiral. We use Eq. (fT5|) 



FIG. 9: (Color online) The AUC, AMD, TFT and WSLS com- 
ponents of the mixed fixed point for the infinitely iterated PD 
for changing u at e = 10~ 4 . Microscopic dynamics are of the 
Moran type. Dashed lines are analytical results from a linear 
expansion in u, circles are from a numerical evaluation of the 
fixed point. The vertical solid black line indicates the loca- 
tion of U W i.e. the value of u where the mixed fixed point 
changes stability. 
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FIG. 10: The stability of the mixed fixed point for Moran 
updating in the infinitely iterated PD game in the (e, u) plane. 
Note that the border between the limit cycle and spiral saddle 
regions is approximate. 



with four strategies to track the location and stability 
properties of the mixed fixed point as u is varied. The 
path of the mixed fixed point at constant e and changing 
u for the Moran process, is shown in Fig. [9] The dashed 
lines are the result of a similar perturbative expansion 
to that carried out in Sec. IIIII where again we see good 
agreement with numerical results for changes in u up to 
u«0.01. 

Similar to our analysis of the three-strategy game, we 
can determine the stability of the mixed fixed point as a 
function of the model parameters e and u. The classifi- 
cation of the nature of the fixed points is more involved 
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FIG. 11: (Color online) The power spectra for the four strate- 
gies with Moran updating at u = 0.02 and e = 0.01. Symbols 
are results from numerical simulations at TV = 12000, solid 
lines are the predictions of Eq. (|31l) . 



for the four-strategy game, however, as we are analyz- 
ing a three-dimensional dynamical system. Stable spirals 
are now fixed points with one pair of complex-conjugate 
eigenvalues with a negative real part and an additional 
real- valued negative eigenvalue. If the sign of the real 
part of the pair of complex-conjugate eigenvalues is op- 
posite to that of the real-valued eigenvalue we will refer 
to the fixed point as a spiral saddle [49]. Fixed points 
with three real- valued negative eigenvalues of the Jaco- 
bian are referred to as stable nodes. The resulting phase 
diagram for the Moran dynamics is shown in Fig.[TUl The 
other three types of microscopic dynamics give qualita- 
tively similar phase diagrams, but the exact quantitative 
positions of the various phase lines will generally be dif- 
ferent. 

When the mixed fixed point is a spiral saddle the deter- 
ministic dynamics can either converge to a limit cycle or 
to the attractor at A11D, depending on initial conditions. 
For locations in the parameter space where the mixed 
fixed point of the deterministic dynamics is a stable spi- 
ral, we again observe demographic oscillations, and they 
can be characterized analytically in a manner similar to 
that discussed in the previous section. We depict the re- 
sulting power spectra for the Moran process in Fig. [TT] 
As seen in the figure WSLS and A11D in particular un- 
dergo strong demographic oscillations, with a comparable 
magnitude between the two strategies. 

As with the iterated PD considered earlier, the am- 
plitude of quasi-cycles resulting from an amplification of 
intrinsic fluctuations can be expected to be large when a 
fixed point of the stable spiral type is located close to the 
border between the stable-spiral and limit-cycle phases. 
There can then again be periods in time when WSLS 
is the most prevalent strategy, despite A11D dominating 
the fixed point. The power spectra for oscillations in the 
concentration of the TFT strategy resulting from the four 
different update rules are compared in Fig. 1121 We again 
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FIG. 12: (Color online) The oscillations in the TFT strategy 
concentration for the four update rules at u — 0.05 and e = 
10 -4 . Symbols are from simulations, lines from the theory. 
System sizes in the range N = 10 4 to 10 6 are used. For 
this particular choice of parameters, the spectra for the local 
and Fermi processes overlap, although this is not the case in 
general. 



observe that the Fermi process exhibits oscillations of a 
higher frequency and with a larger amplitude than the 
other rules. Although spectra for the local process and 
the Fermi process overlap in the figure, this is coinciden- 
tal at this point in parameter space, and will not be the 
case in general. 



VI. CONCLUSION AND OUTLOOK 

In this paper we have used analytical approaches based 
on van Kampen's system-size expansion to study the 
emergence of quasi-cycles in evolutionary games in finite 
populations. Most existing studies of such effects are of 
a numerical nature; we have complemented these giving 
a systematic account of the formalism, and derived the 
resulting effective Langevin equations which describe the 
statistics and correlations of fluctuations in large, but fi- 
nite populations. This approach and our general formu- 
lae are applicable to a large class of microscopic update 
rules, and to games with an arbitrary number of pure 
strategies. They are, in principle, also valid for arbitrary 
mutation matrices. The results of this paper hence allow 
one to predict the regions in parameter space in which 
coherent quasi-cycles are to be expected, and to compute 
their spectral properties. In particular coherent cycles, 
such as reported in game dynamical systems e.g. in [l4j . 
can be understood as a consequence of the combination 
of intrinsic noise and the existence of a stable fixed point 
with complex eigenvalues in the corresponding determin- 
istic system obtained in the limit of infinite populations. 
In absence of noise the deterministic system approaches 
such fixed points in an oscillatory manner, with oscilla- 
tions dying away at a rate proportional to the inverse 
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real part of the relevant eigenvalue. In finite systems, 
discretization noise leads to persistent stochastic correc- 
tions perturbing the system at all times. In the limit of 
large, but finite system sizes these fluctuations (the noise 
x) in Eq. (f2T)|). together with the pre-factor TV -1 / 2 ) can be 
seen as a small perturbation to the deterministic dynam- 
ics, driving the system away from the fixed point. The 
attracting fixed point and the oscillatory approach to it 
on the deterministic level on the one hand and the per- 
sistent intrinsic noise on the other then conspire to give 
coherent and sustained stochastic cycles, with an ampli- 
tude largely determined by the inverse real part of the 
least stable complex eigenvalue. 

We have applied the van Kampen formalism to the spe- 
cific example of the iterated PD game, where stochastic 
oscillations have been reported in the earlier numerical 
study [HI • We have worked out detailed phase diagrams 
depicting the nature of the limiting deterministic dynam- 
ics and we have studied systematically how the muta- 
tion rate and complexity cost, the two main model pa- 
rameters, affect the outcome of the deterministic system. 
Based on this analysis we are able to predict the parame- 
ter regimes in which stochastic oscillations occur. In par- 
ticular we find that oscillation amplitudes become maxi- 
mal when the Hopf bifurcation line in the phase diagram 
is approached from within the stable phase. At the bi- 
furcation line the oscillation amplitude formally diverges, 
with the power spectrum turning into a delta-function, 
and as the instability line is crossed a limit cycle is born. 

We have also carried out a detailed comparison of four 
different microscopic update rules; results indicate that 
their respective phase diagrams are qualitatively similar. 
The analysis shows that, at fixed values of the model pa- 
rameters, the Fermi process tends to produce stochastic 
cycles with larger amplitudes and frequencies than the 
other update rules. We have extended our study to a 
version of the iterated PD game in which errors of the 
trembling-hand type occur with a small, but non-zero 
rate. The so-called win-stay lose-shift strategy has here 
been seen to out compete tit-for-tat, and accordingly we 



have considered a four-strategy space (always defect, al- 
ways co-operate, tit-for-tat and win-stay lose-shift), and 
have identified the regions in parameter space where co- 
herent cycles are most likely to occur. Analytical results 
for the resulting power spectra of these quasi-cycles are 
confirmed convincingly in numerical simulations. 

Mathematical techniques of the type we have used 
here, most notably the master equation formalism and 
system-size expansions, were first devised in statistical 
physics, but they are becoming increasingly more popu- 
lar in the game theory literature. This extends to equiv- 
alent approaches based on Kramers-Moyal expansions. 
We attribute this popularity to the generality with which 
these methods are applicable and to the fact that they 
allow one to obtain an exhaustive account of the prop- 
erties of first-order stochastic corrections to the limiting 
deterministic dynamics. Exact analytical results can be 
derived for large, but finite populations, and hence these 
techniques make simulations on the microscopic level re- 
dundant (at least in principle). We expect this to be 
the case for games with more complicated strategy struc- 
tures, or with interaction between more than two players 
such as for example in public goods games. The analyt- 
ical approach can also be expected to be applicable to 
other, potentially more intricate types of human error. 
For such games it may be difficult or time consuming 
to carry out reliable simulations and to perform exhaus- 
tive parameter searches. Analytical characterizations of 
stochastic effects such as the ones discussed in this paper 
may then be particularly welcome. 
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Appendix A: System-size expansion of the master 
equation 

The van Kampen system-size expansion has been ex- 
tensively discussed elsewhere, together with explicit ex- 



amples [6H9> and so here we will only briefly summarize 
the general idea and give the final results of the calcula- 
tion that are relevant for this paper. 

The starting point is the substitution of the ansatz ([27)1 
into the master equation — or its generalization to 
more than three strategies. This yields an expansion in 
powers of N^ 1 / 2 , after a re-scaling of time by a factor 
of N . To leading order (N~ x / 2 ) the deterministic equa- 
tion (fl3|) is obtained. To next-to-leading order (iV _1 ) 
the Fokker-Planck equation ([25)1 is found. This is de- 
fined in terms of two quantities: C,; = J^k Jik£,k, where 
J is the Jacobian of the dynamics given by Eq. (fl~3|) . and 
B a symmetric matrix. Since we are interested in fluc- 
tuations about stationary states, both and -B,j are 
time- independent . 

The Jacobian can be obtained in a straightforward 
fashion once the dynamics (fTBf is known. The elements 
of the matrix B are found from the iV _1 terms in the 
system-size expansion to be 




-[T i ^ j {x)+T j ^ i (x)] , iii^j 



_ . . (AD 

Therefore the deterministic and stochastic dynamics to 
the order we are working are entirely determined by 
Ti^j(x). This can be found by making the substitu- 
tions (m/N) -> x u fi -> f°° and (f) -> 0°° in Eq. (fT0|) . 
to obtain 

lU^x) = ^x k x i g ki {f°°)q kj . (A2) 

k 

This explicitly shows how to construct (x) , once the 
process (defined by g k i{f)) and the mutation matrix (q%j) 
have been given. 



Appendix B: Payoff matrix for a four strategy, 
infinitely repeated PD with trembling hand errors 

When two players meet and play the PD game over 
multiple rounds their state in round I is defined by their 
actions in that round, e.g. player 1 co-operates, player 
2 defects. The actions of the pair of players then deter- 
mines the payoff they each receive. The payoff matrix 
for the infinitely repeated PD game with trembling hand 
errors is constructed by considering the stationary distri- 
butions of the state of each pair of players, to first order 
in e, the probability of a 'trembling hand' error occur- 
ring [23} . It is given by 
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AUC AUD TFT WSLS 



AUG ( R-e(2R- S-T) S + e(R + P-2S) R — e(3P — T — 2S) (R + S)/2 + (e/2)a 

AUD T-e(2T-R-P) P + e(S + T-2P) P + e(S + 2T - 3P) (P + T)/2 - (e/2)a 

TFT R + e(2T + S-3R) P + e(T + 2S - 3P) 7 7 

WSLS 1 \ (P + T)/2 + {e/2)(3 (P + S)/2 - (e/2)f3 7 P + e(T + 2P + - 4P) 




where a = (T + P - R - S), /3 = (S + P - R - T) and 7 = (r + P + P + 5)/4. 



